function [first_spike, first_isi, nr_spikes, spike_train] = find_spikes(DFILE, tsw1i, tsw2i, spike_threshf)
%finds the latency of the first spike from t_zero
% voltage is the voltage matrix (records, voltage) from datac2mat
% rates can change and are reported by DFILE
%
% tsw1 and tsw2 are the times for the spike windows begin and end
%
% Modified 09 Feb 2000 P. Manis
%  - take rate array (each record, sync'd to voltage matrix) and times to get window for analysis
%    in case sampling rate changes in middle of analysis (why not?)
%
% Changed to pass DFILE - the only elements of DFILE that are accessed here are
% DFILE.rate and DFILE.nr_channel. (so you could easily fake it).
% 
% Modified 4/19/2002 P. Manis
% slightly faster detection algorithm - also uses less memory.
%
global VOLTAGE
intermediate_plot = 0;

[recs,pts]=size(VOLTAGE);

if(length(spike_threshf) > 1) % handle events in both directions....
   spike_thresh = spike_threshf(1);
   sign = -1;
else
   spike_thresh = spike_threshf;
   sign = 1;
end;

   
% initialize the output arrays to null - default is we find no spikes...
first_spike = [];
first_isi = [];
nr_spikes = 0;
spike_train = [];
sthr = 1;
k = 0;
RATES = DFILE.rate.*DFILE.nr_channel/1000;
time = make_time(DFILE);
if(length(tsw1i) == 1)
   tsw1 = tsw1i*ones(1, recs);
   tsw2 = tsw2i*ones(1, recs);
else
   tsw1 = tsw1i;
   tsw2 = tsw2i;
end;


for i=1:recs
	sw1 = floor(tsw1(i)/RATES(i));
   if(sw1 < 1)
      sw1 = 1;
   end;
   if(sw1 > length(VOLTAGE))
      sw1 = length(VOLTAGE-1);
   end;
   sw2 = floor(tsw2(i)/RATES(i));
	if(sw2 > length(VOLTAGE))
		sw2 = length(VOLTAGE)-1;
	end;
% find spikes
   k1=zeros(1, length(VOLTAGE(i, sw1:sw2))); % initialize pulse vector
   if(sign > 0)
      [d1, k2] = find(VOLTAGE(i,sw1:sw2)>spike_thresh); % find points > threshold
      switch(length(k2))
      case 0
         spike_pos = [];
      case 1
         spike_pos = k2(1);
      otherwise
         spike_pos = [k2(1) k2(find(diff(k2) > 1)+1)]; % and store the resulting array
      end;
      
   else
      [d1, k2] = find(VOLTAGE(i,sw1:sw2)<spike_thresh); % find points > threshold
      switch(length(k2))
      case 0
         spike_pos = [];
      case 1
         spike_pos = k2(1);
      otherwise
         spike_pos = [k2(1) k2(find(diff(k2) > 1)+1)]; % and store the resulting array
      end;
   end;
   % now transform spike latency from sample time to actual latencies
   % this involves interpolating the trace.
   spike=[];
   for j=1: length(spike_pos)
       vx = VOLTAGE(i, sw1+spike_pos(j)-1);
       m = vx - VOLTAGE(i, sw1+spike_pos(j)-2);
       b = vx - m*spike_pos(j); % compute the intercept (note, using points at time units...)
       if(m == 0) % points were the same... just take one
           spike(j)=spike_pos(j);
       else
           spike(j) = ((spike_thresh - b) / m); % get corresponding point in interval
       end;
   end
   spike = spike * RATES(i);
   if(i == recs && intermediate_plot)
      figure
      ld1 = length(d1);
      ld2=length(d2);
      ld3=length(d3);
      hold on;
      %        plot(time(1:ld1), d1, 'r');
      plot(time(i, 1:ld2), d2, 'g');
      plot(time(i, 1:ld3), d3, 'b');
      hold off;        
   end
   nr_spikes(i)=length(spike);
   if nr_spikes(i) >= 1
      first_spike(i)=spike(1);
   else
      first_spike(i)=NaN;
   end
   if nr_spikes(i) > 1
      first_isi(i)=spike(2)-spike(1);
   else
      first_isi(i)=NaN;
   end
   if(nr_spikes(i) > 0)
      k=k+1;
      spike_train(k).latency = spike; % save each spike train for the latency/isi plot
      spike_train(k).source = i; % save the source index (which voltage trace) with the spike 
   end
end




